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Abstract. After reviewing the main features of anomalous energy transport in 
ID systems, we report simulations performed with chains of noisy anharmonic 
oscillators. The stochastic terms are added in such a way to conserve total energy 
and momentum, thus keeping the basic hydrodynamic features of these models. 
The addition of this " conservative noise" allows to obtain a more efficient estimate 
of the power-law divergence of heat conductivity k{L) ~ L°' in the limit of small 
noise and large system size L. By comparing the numerical results with rigorous 
predictions obtained for the harmonic chain, we show how finite-size and -time 
effects can be effectively controlled. For low noise amplitudes, the a values are 
close to 1/3 for asymmetric potentials and to 0.4 for symmetric ones. These results 
support the previously conjectured two-universality-classes scenario. 

1 Introduction 

Many-body systems constrained in reduced spatial dimensions (1 and 2D) display unusual 
relaxation and transport properties. Following familiar arguments of equilibrium statistical 
mechanics, this should be traced back to the predominant role of fluctuations that give rise 
to, e.g., long-range order etc. In the context of non-equilibrium processes, the presence of 
anomalies is signalled by the appearance of long-time tails in the correlation functions of 
the relevant currents 1 2 leading to ill-defined transport coefficients i.e. to the breakdown of 
standard hydrodynamics. The case of one-dimensional models is perhaps the most striking. 
For the sake of an example, we mention the divergence of viscosity in cellular automata fluids 
[3], anomalous diffusion in single-file systems 4J, and the enhancement of vibrational energy 
transmission in polymers or individual carbon nanotubes |^. 

Within this general context, one of the issues that attracted a renewed interest in the last 
decade is the problem of anomalous heat conduction in one-dimensional models (for a review 
see [7]). The interest in such models originates from the need of constructing a minimal, nonper- 
turbative theory of nonequilibrium stationary states and the quest for a rigorous microscopic 
foundation of phenomenological relations (the Fourier's law in this case). This motivation gen- 
erated a vast literature, expecially in the mathematical physics community [S]. Moreover, many 
of the peculiarities of ID models turned out to be of interest by themselves, as examples of 
highly correlated, and thus complex, behaviour. Those latter features are shortly reviewed in 
the first part of this paper. In the second part, the usefulness of an additional stochastic noise 
[9] is discussed with reference to some open questions. 
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The class of models we will consider consists of a set of N classical point -like particles with 
masses r7i„ and positions x„ ordered along a line. Interactions are restricted only to nearest- 
neighbour pairs and the dynamics is ruled by 

rrinXn = -Fn + F„_i ; Fn = -V'{Xn+l - Xn) , (1) 

where V'{y) is a shorthand notation for the first derivative of the the interparticle potential 
V with respect to the argument. Boundary conditions for the first and last particle will be 
specified in the various cases. The two most intensively studied systems are the Fermi-Pasta- 
Ulam (FPU) model (with equal masses to„ = m) |10|ll|12|13j 0, 

V{y) = 'f{y-ar + f{y-a)^ + '-±{y-a)\ (2) 

and the diatomic Hard- Point Gas (HPG), where m„ — m (rm) for even (odd) n, with the 
interaction potential [14115116117] . 

~ 1 otherwise 

We clarify from the beginning that we will always deal with genuine nonintegrable dynamics. 
For the FPU, this means working with high enough energies/temperatures to avoid all the 
difficulties induced by quasi-integrability and the associated slow relaxation to equilibrium 
jT8j . For the HPG this requires fixing a mass-ratio r not too close to unity. 

The results that emerged from a long series of works can be summarized as follows. All 
studied models of the form ([1]) display anomalous transport and relaxation features, meaning 
by this that (at least) one of the following phenomena has been reported: 

— The finite-size heat conductivity k{L) diverges as in the limit of a large system size 
L —> oo |13) . This means that this transport coefficient is ill-defined in the thermodynamic 
limit, i.e. Fourier's law does not hold; 

— The equilibrium correlator of the energy current displays a nonintegrable power-law decay, 
(J(t)J(O)) cx (0 < (5 < 1) for long times t ^ oo [TOj. Accordingly, the Green-Kubo 
formula yields an infinite value of the conductivity; 

— Energy perturbations progate superdiffusively [20j : a local perturbation of the energy broad- 
ens and its variance cr^(i) grows in time as t^ with /3 > 1; 

— Relaxation of spontaneous fluctuations is fast (i.e. superexponential) |21I22| : at variance 
with standard hydrodynamics, the typical decay rate in time of fluctuations, T(q), is found 
to scale as r(q) (with z < 2). 

Altogether, these features indicate that the kinetics of energy carriers is so correlated that 
they are able to propagate faster with respect to the the standard (diffusive) case. In view 
of this common physical origin, it is expected that the exponents describing such process will 
be related to each other by some "hyperscaling relations" . Their value should be ultimately 
dictated by the dynamical scaling of the underlying dynamics. Moreover, we may at least hope 
that they are largely independent of the microscopic details, thus allowing for a classification 
of anomalous behaviour in terms of "universality classes" . 

Numerical studies [7] indicate that anomalies occur generically in 1 and 2D, whenever 
momentum is conserved. This is connected to the existence of long-wavelength (Goldstone) 
modes (an acoustic phonon branch in the linear spectrum of ([1])) that are very weakly damped. 
Indeed, it is sufficient to add external (e.g. substrate) forces to make the anomaly disappear. 
This is precisely the case of the ding-a-ling model [23] and of other models in the same class, like 
the Frenkel-Kontorova |24I25| or the nonlinear Klein-Gordon chains [55]. The only remarkable 
exception is the coupled rotor chain [27|28| where, however, different mechanisms are at work. 

^ For the sake of nomenclature, it is useful to mention two particular cases: the quadratic plus cubic 
((74 ~ 0) and quadratic plus quartic {gs — 0) potentials that, for historical reasons, are referred to as 
the "a-FPU" and "/3-FPU" models, respectively. In the former one, sufficiently small coupling constant 
(/3 and/or energies must be considered to avoid runaway instabilities. 
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Recently introduced stochastic models provide some mathematically rigorous results about 
the importance of momentum conservation. A random exchange of momentum between neigh- 
bor particles is added to the Hamitonian deterministic dynamics. These exchanges may conserve 
only the energy [29|, or also the total momentum [9]. While in the first case conductivity re- 
mains finite and Fourier law is proved for harmonic interaction [29j . when total momentum 
is conserved it diverges as L^/^ [9 . In the anharmonic case, conductivity is much harder to 
compute or estimate also with this noise. So we performed some numerical simulations on the 
quartic and cubic FPU models with momentum-|-energy conserving noise. 

In the first part of the paper we will review in more detail the anomalous features mentioned 
above. The second part is devoted to some preliminary numerical studies of the stochastic 
models. 



1.1 Diverging finite-size conductivity 

A natural way to simulate a heat conduction experiment consists in putting the system in con- 
tact with two heat reservoirs operating at different temperatures r+ and r_. Several models 
for the thermostats have been proposed based on both deterministic and stochastic algorithms 
[?]• Regardless of the actual thermostatting scheme, after a transient, an off-equilibrium sta- 
tionary state sets in, with a net heat current flowing through. The thermal conductivity k 
is then estimated as the ratio between the time-averaged flux j and the overall temperature 
gradient (r+ — T^)/L, where L — aN is the chain length (a denoting the lattice spacing). In 
this manner, k should be considered as an effective transport coefficient, accounting for both 
boundary and bulk scattering mechanisms. The average j can be estimated in several equivalent 
ways, depending on the employed thermostatting scheme. One possibility is to directly measure 
the energy exchanges with the two baths. A more general definition (thermostat-independent) 
consists in averaging 

jn ^ -^{in+l + in) F{Xn+l ~ Xn) , (3) 

which is obtained from the continuity equation for the energy density 



_ 1 .2 1 

— Jn — 1 Jn 



V{Xn+l - Xn) + V{Xn - Xn-l) 



(4) 



As a result of many independent simulations performed with the above-described methods, 
it is now established that k (x N". It is also remarkable that the same type of behaviour 
has been observed for a realistic model of a single-walled carbon nanotube 6 . This type of 
molecular dynamics simulations, that involves complicated three-body interactions for carbon 
atoms, confirm that toy models like ours can indeed capture some general features. 



1.2 Long-time tails 

In the spirit of linear-response theory, transport coefficients can be computed from equilibrium 
fluctuations of the associated currents. More precisely, by introducing the total heat flux 

^ = EJ"' (5) 

n 

the Green-Kubo formalism tells us that heat conductivity is determined from the expression 

^ lim lim 4 / dt' {J{t')J{0)) (6) 



where the average is performed in a suitable equilibrium ensemble, e.g. microcanonical with 
zero total momentum. 
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A condition for the formula ([6]) to give a well-defined conductivity is that the time integral 
is convergent. This is clearly not the case when the current correlator vanishes as ( J(t) J(0)) (x 
^-(1-^) with < S < 1. Here, the integral diverges as and we may thus define a finite-size 
conductivity k{L) by truncating the time integral in the above equation to i « L/c, where c 
is the sound velocity. Consistency with the above definition thus implies a = S. The available 
data are consistent with this expectation, thus providing independent method for estimating 
the exponent a. 

For later purposes, we mention that, by means of the Wiener-Khintchine theorem, one can 
equivalently extract S from the low-frequency behaviour of the spectrum of current fluctuations 

Siu;) = j dij(J{t)J{0))e"^' (7) 

that displays a low-frequency singularity of the form S{uj) cx uj^^ . From the practical point of 
view, this turns out to be the most accurate numerical strategy as divergence rates are better 
estimated than vanishing ones. 



1.3 DifFusion of perturbations 

Consider the infinite system at equilibrium with an energy eo per particle and average mo- 
mentum 0, and perturb it by increasing the energy of a subset of adjacent particles by some 
preassigned amount Ae . Let us denote by e{x,t) the energy profile evolving from such a per- 
turbed initial condition (for simplicity we identify x with the average particle location ni) . We 
then ask how the perturbation 

f{x,t)^ {e{x,t)-eo) (8) 

behaves in time and space j30| . where the angular brackets denote an ensemble average over 
independent trajectories. Because of energy conservation, J2n fi''^^^^) ~ times, so 

that / can be interpreted as a probability density (provided that it is also positive-defined and 
normalized) . 

At sufficiently long times and for large x, one expects f{x,t) to scale as 

/(x,t)=i-^^(xA^) (9) 

for some probability distribution Q and a parameter < 7 < 1- 

The case 7 = 1/2 corresponds to normal diffusion and to a normal conductivity. On the other 
hand, 7=1 corresponds to a ballistic motion and to a linear divergence of the conductivity. 

Consequently an a strictly contained between and 1 implies a superdiffusive behaviour of 
the macroscopic evolution of an energy perturbation. It is an open problem to determine the 
dynamical nature of this superdiffusion. In [31j has been proposed that Levy walks 132] may 
describe these dynamics. 



1.4 Relaxation of spontaneous fluctuations 

The evolution of a fluctuation of wavcnumber q excited at t = is described by its corre- 
lation functions G{q, t) . For ID models like llj they are deflned by considering the relative 
displacements m„ = a;„ — and defining the collective coordinates 

1 2-Kk 
QiQ) = -7^ ^" exp(-ign) 9 = ^ (10) 

^ n—l 

for k being an integer comprised between —N/2 + 1 and N/2. The normalized correlator 
G{q,t) = ((5*((7, i)(5(g, 0))/(|(5(g')p), is thus proportional to the density-density correlator, 
which is routinely used in condensed-matter physics |33j . 
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For sufficiently small wavenumbers q, the Fourier transform of G, the dynamical structure 
factor S{q,uj), usually displays sharp peaks at finite frequency. The associated linewidths are 
a measure of the fluctuation's inverse lifetime. Simulations 18 22] indicate that these lifetimes 
scale as q~^ with z = 1.50 — 1.67. As explained above, one may think of this as a superdiffusive 
process, intermediate between standard diffusive and ballistic propagation. 

This property is supported by theoretical results obtained in the framework of Mode- 
Coupling Theory (MCT). This approach has been traditionally invoked to estimate long-time 
tails of fluids [i^ and to describe the glass transition [36,^ . Basically, it amounts to writing a set 
of approximate equations for G{q, t) that must be solved self-consistently. For the problem at 
hand, the simplest version of the theory amounts to consider the equations [37|38j 

G{q,t) + e f r{q,t^ s)G{q,s)ds + uj\q)G{q,t) = , 
Jo 

where the memory kernel r{q,t) is proportional to {J-{q,t)T{q,0)), with being the non- 
linear part of the fluctuating force between particles. Equations (fTTjl must be solved with the 
initial conditions G{q, 0) = 1 and G(q, 0) = 0. Equations (fTT|) are derived within the well-known 
Mori-Zwanzig projection approach [34j . 

The mode-coupling approach basically amounts to replacing the exact memory function F 
with an approximate one, where higher-orders correlators are written in terms of G{q, t). In the 
generic case, in which 173 is different from zero, the lowest-order mode coupling approximation 
of the memory kernel turns out to be |37|38j 

r{q,t)^^^{q)^-^ Gip,t)Gip\t) . (11) 

Here p and p' range over the whole Brillouin zone (from — tt to tt in our units) . This yields 
a closed system of nonlinear integro-differential equations. Both the coupling constant e and 
the frequency uj{q) are temperature-dependent input parameters, which should be computed 
independently by numerical simulations or approximate analytical estimates [37|38) . For the 
aims of the present work, we may restrict ourselves to considering their bare values, obtained 
in the harmonic approximation. In the adopted dimensionless units they read e = 3f/|fcBT/27r 
and uj{q) = 2|sin||. Of course, the actual renormalized values are needed for a quantitative 
comparison with specific models. 

The long-time behaviour of G can be determined by looking for a solution of the form 

G{q,t) = C(g,t)e^"(')*+c.c. . (12) 

with G <C ojG. It can thus be shown |39|40j that, for small q- values and long times C{q,t) ~ 
g{^/£tq^^^) i.e. z — 3/2 in agreement with the above mentioned numerics. Furthermore, in the 
limit ^tq^/"^ one can explicitely evaluate the functional form of g, thus obtaining 

C{q,t) = ^exp{-Dq'\t\^) , (13) 

where Z? is a suitable constant of order unity. Correlation display a "compressed exponential" 
behaviour [IT] in this time range. This also means that the lineshapes of the structure factors 
S{q,u)) are non-Lorenzian but rather display a faster power-law decay (w — uJmax)~'^^^ around 
their maximum. 

Upon inserting this scaling result into the definition of the heat flux, one eventually concludes 
that the conductivity k diverges with a rate a — 1/3. 

2 Universality 

The crucial question at this point is: how universal are the above defined exponents? The 
renormalization group argument by Narayan and Ramaswamy [42], predicts a = {2 — d)/{2 + d). 
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Following the arguments exposed above, this implies that in ID the values of the exponents are 

« = ^ = ^, z = = ^ . (14) 

According to this approach, any possible additional term in the noisy Navier-Stokes equation 
yields irrelevant corrections in the renormalization procedure j42j . 

On the other hand, there now exists a rigorous result, proving that a = 1/2 [9J in a chain 
of harmonic oscillators subject to an energy and momentum-conserving noise. Moreover, the 
application of kinetic theories to the /3-FPU model ^43,44, 45] predict a — 2/5, while a modified 
version of the MCT adapted to this specific case gives instead a = 1/2 [46l40j . 

The validation of these theoretical results by numerical simulations is, to some extent, 
challenging. Generally speaking, the available numerical estimates of a and S range between 
0.25 and 0.44 [7148] . The existence of crossovers among different scaling regimes has been 
observed [50j . However, even in the most favorable cases of computationally efficient models as 
the HPG, finite-size corrections to scaling are sizeable. As a matter of fact, a-values as diverse 
as 0.33 [17] and 0.25 [49^ for comparable parameter choices have been reported. On the other 
hand, a numerically convincing confirmation of the a = 1/3 prediction comes from the diffusion 
of perturbations [2U] . The most compelling deviations from the values given in (fH)) have been 
reported for the /3-FPU model [48] , where a better agreement with the predictions of the kinetic 
theory has been observed. 

A reasonable argument that can be invoked to delimit the a = 1/3 universality class appears 
to be the symmetry of the interaction potential with respect to the equilibrium position. In 
fact, systematically larger a-values have been reported only for symmetric potentials such as 
the /3-FPU model or, more recently, for a modified version of the HPG model [3S], where the 
interaction can be tuned to yield symmetric fluctuations of the force. Moreover, the rigorous 
result a = 1/2 [9] indeed concerns a model whose potential (harmonic) is symmetric and where 
the stochastic updating is symmetric as well. In the framework of MCT, one can understand 
that the symmetry of the fluctuations implies that the quadratic kernel ([TT|) should be replaced 
by a cubic one [?(T , thus yielding different values of the exponents (see [33] for a thermodynamic 
interpretation of this difference). 

Whether symmetry is the only necessary ingredient to identify the asymptotic behaviour is 
however not fully proven. As mentioned above, difficulties manifestly arise e.g. for the /?-FPU 
model where it is not yet completely assessed which of the predicted exponents one should 
expect. In this respect some further reconsideration of MCT and kinetic theories may be of 
help. Actually, there are even some controversial results about the dynamical scaling exponent 
z in the supposed broad a = 1/3 class [47]. It is thus clear that further rigorous results on 
simple models would be definitely welcome. 



3 Stochastic models with energy and momentum conservation 

Recently, some of the authors proposed a new class of models for analyzing the anomalous 
properties of heat conduction for a system of oscillators [S]. The Hamiltonian dynamics is 
perturbed by a stochastic noise, which acts only on the momenta, while conserving total energy 
and momentum. In particular, such an approach allows to compute explicitly the heat flux 
correlation function in the harmonic case, i.e. model ^ with 33 = 54 = 0. For the ID model it 
is found that the current correlator vanishes as (J{t)J{0)) ^ t^^^^, thus implying that k ^ L^^^. 
The average of the current correlator is performed over the equilibrium measure, which is the 
uniform measure over the hypersurface of constant total energy and momentum. As usual 
in this microcanonical measure, the total momentum is set to zero. There are various ways 
for translating such stochastic dynamics into a suitable algorithm. The numerical simulations 
hereafter reported have been performed by making the system evolve over a finite number n 
of steps (each of duration h) and, then, by updating the momenta of triplets of nearest- 
neighbor oscillators, whose position are randomly chosen with a uniform probability density 
over the chain sites. Each "collision" in the momentum-subspace is performed by extracting 



Will be inserted by the editor 



7 



from the uniform distribution over [0, 27r] a random angle which rotates the three momenta in 
such a way to maintain fixed their sum and the sum of their squares. The resulting configuration 
is the initial condition for a new deterministic trajectory lasting again over n steps, and so forth. 
Because of the conservation laws, the microcanonical measure (uniform on the energy surface) 
is still invariant and it can be proven to be ergodic. 

The algorithmic procedure introduces a new control parameter in the dynamics, namely the 
ratio between the number of collisions per the unit time and the number of particles N . 

In order to test the effectiveness of the algorithm, we have first performed numerical sim- 
ulations of the stochastic dynamics of the harmonic chain, as we can make a comparison with 
the rigorous results . In the left panel of figure [T] we show the power spectrum defined in ([7]) 
for different values of N , while maintaining fixed the ratio ric/N ~ 0.1. All curves align to the 
same power-law behavior in a range of small values of uj of about 3 decades. In particular, the 
spectrum is found to diverge as S{uj) ^ uj~^ with S slightly larger than 1/2. Moreover, finite size 
effects are clearly visible in the spectrum corresponding to the smallest chain length {N = 512). 
In the right panel of figure [l] we have kept N = 2048, while we have considered different values 
of ric- For decreasing values of ric, we see that 6 remains close to 1/2, in agreement with the 
rigorous prediction reported in [9] . 




10'^ lO''* 10'' 10'^ 10'"' 10° 1^^ 10"^ lO"'' 10"' 10"^ 10"' 10° 1^' 



Fig. 1. Power spectra of the heat flux for the stochastic harmonic chain at energy density e = 10. Units 
along the vertical axis are arbitrary. Left panel: data for different values of A'^, while maintaining fixed 
the ratio ric/N. Right panel: A*' = 2048, and different values of ric. Each curve is averaged over about 
500 independent trajectories. Here and in the following figures, the time interval between collisions is 
10 h with integration time step h = 0.01. The thin dashed line corresponds to the law w"^'''^. 



Relying on these results, we have also studied the power spectra of the stochastic FPU 
models ([2]) with g2 ~ gs = 54 = 1 and g2 — ~ 1 and g^ = (/3-FPU). For both models, 
we report the power spectrum at fixed N for different Uc values (see Figs. I2I3|) . Analogously 
to what observed for the harmonic case, we see that upon increasing ric, the scaling region 
widens towards higher frequencies, thus allowing for a more accurate analysis than in the 
purely deterministic case. Moreover, in both cases the effective divergence exponent S appears 
to increase with ric, approaching 1/2. More specifically, for the FPU with cubic and quartic 
nonlinearities (Fig. [2]), S increases from « 0.35 for ric = 10 to « 0.48 for Uc = 200. Analogously, 
for the /3-FPU case (Fig. [3]), we find that 3 increases from « 0.41, to « 0.47. 

Based on the existence of a few universality classes, a continuous dependence of the exponent 
6 on the number of collisions is quite unlikely. However, one should recognize that the 
frequency scaling range is as wide as three decades and the dependence oi S on uj is pretty 
weak. Altogether, the fact that S is close to 1/2 for the larger Uc values suggests that the scaling 
behaviour predicted in [S] applies to a larger class of systems and enforces the hypothesis of a 
second distinct universality class. Moreover, it is interesting to notice that in both anharmonic 
models, the divergence of thermal conductivity increases with the strength of the noise. This 
means that, quite surprisingly, noise contributes to slowing down the decay of energy current 
correlations. 
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Fig. 2. Power spectrum of the heat flux for the stochastic FPU model (cubic and quartic nonlinearities) 
with energy density e = 10, = 2048 and different number Uc of colliding triplets. Units along the 
vertical axis are arbitrary. Each curve is averaged over about 800 independent trajectories. The thin 
dashed line corresponds to the law 




-4 -? -1 n 1 
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Fig. 3. Power spectrum of the heat flux for the stochastic /3-FPU model (quartic nonlinearity) with 
energy density e = 10, A'^ = 2048 and different number Uc of colliding triplets. Units along the vertical 
axis are arbitrary. Each curve is averaged over about 800 independent trajectories. The thin dashed 
line corresponds to the law oj"^''^. 
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